This creates points at the geometric center of each polygon.
Code
mapview(Centroid_geo)
Code
mapview(Centroid_geo) +mapview(Municipalities_geo, alpha.regions =0) # both maps, with full transparency in polygons
But… is this the best way to represent the center of a transport zone?
These results may be biased by the shape of the polygons, and not represent where activities are. Example: lakes, forests, etc.
To overcome this, we can use weighted centroids.
10.2 Weighted centroids
We will weight centroids of transport zones by population and by number of buildings.
For this, we will need the census data (INE 2022).
Code
Census =st_read("../data/census.gpkg", quiet =TRUE)mapview(Census |>filter(Municipality =="Lisboa"), zcol ="Population")
It was not that easy to estimate weighted centroids with R, as it is with GIS software. But there is this new package centr that can help us (Zomorrodi 2024).
We need to specify the group we want to calculate the mean centroids, and the weight variable we want to use.
Code
library(centr)Centroid_pop = Census |>mean_center(group ="Municipality", weight ="Population")
We can do the same for buildings.
Code
Centroid_build = Census |>mean_center(group ="Municipality", weight ="Buildings")
Now, to use plot() we incrementally add layers to the plot.
Code
# Plot the Municipalities_geo polygons first (with no fill)plot(st_geometry(Municipalities_geo), col =NA, border ="black")# Add the Centroids_geo points in blueplot(st_geometry(Centroid_geo), col ="blue", pch =16, add =TRUE) # add!# Add the Centroid_pop points in redplot(st_geometry(Centroid_pop), col ="red", pch =16, add =TRUE)# Add the Centroid_build points in blackplot(st_geometry(Centroid_build), col ="black", pch =16, add =TRUE)
Static map of different centroids of Municipalities
In the next section we will use these centroids to calculate the desire lines between them.
MQATMQAT11OD pairs and desire lines9Introduction to spatial data1Introduction2SoftwareData and Models3R basics4Data manipulation5Exploratory Data Analysis6Multiple Linear Regression7Factor Analysis8Cluster AnalysisGIS in R9Introduction to spatial data10Centroids of transport zones11OD pairs and desire lines12Euclidean and routing distances13Open transportation data14Grids15Interactive mapsReferencesGIS in R10Centroids of transport zones
R Félix, G Valença, F Moura
MQAT - Instituto Superior Técnico
September 2025
10Centroids of transport zones – MQAT10Centroids of transport zones – MQAT10Centroids of transport zones – MQATMQAT
---format: pdf: prefer-html: true---# Centroids of transport zonesIn this section we will calculate the geometric and the weighted centroids of transport zones.## Geometric centroidsTaking the `Municipalities_geo` data from the previous section, we will calculate the geometric centroids, using the `st_centroid()` function.```{r geocentroid}#| message: false#| warning: falselibrary(dplyr)library(sf)library(mapview)Municipalities_geo = st_read("../data/Municipalities_geo.gpkg", quiet = TRUE)Centroid_geo = st_centroid(Municipalities_geo)```This creates points at the geometric center of each polygon.quarto-executable-code-5450563D```r#| message: false#| fig-format: pngmapview(Centroid_geo)mapview(Centroid_geo) +mapview(Municipalities_geo, alpha.regions =0) # both maps, with full transparency in polygons```But...is this the best way to represent the center of a transport zone?These results may be biased by the shape of the polygons, and not represent where activities are.Example: lakes, forests, etc.To overcome this, we can use **weighted centroids**.## Weighted centroidsWe will weight centroids of transport zones by population and by number of buildings.For this, we will need the census data [@INEcensus].quarto-executable-code-5450563D```r#| eval: false#| include: false# data prepCENSUS =st_read("/media/rosa/Dados/GIS/MQAT/original/BGRI21_LISBOA.gpkg")CENSUS = CENSUS |>select(DTMN21, SUBSECCAO, N_INDIVIDUOS, N_EDIFICIOS_CLASSICOS) |>rename(Mun_code ="DTMN21",UID ="SUBSECCAO",Population ="N_INDIVIDUOS",Buildings ="N_EDIFICIOS_CLASSICOS")CENSUS = CENSUS |>left_join(Municipalities_names |>select(Mun_code, Municipality) |>distinct()) |>select(Municipality, UID, Population, Buildings)st_write(CENSUS, "../data/census.gpkg", delete_dsn =TRUE)``````{r getcensus}#| message: false#| fig-format: pngCensus = st_read("../data/census.gpkg", quiet = TRUE)mapview(Census |> filter(Municipality == "Lisboa"), zcol = "Population")```It was not that easy to estimate weighted centroids with R, as it is with GIS software.But there is this new package [`centr`](https://ryanzomorrodi.github.io/centr) that can help us [@centr].We need to specify the **group** we want to calculate the **mean centroids**, and the **weight variable** we want to use.```{r centr}library(centr)Centroid_pop = Census |> mean_center(group = "Municipality", weight = "Population")```We can do the same for buildings.```{r buildings}Centroid_build = Census |> mean_center(group = "Municipality", weight = "Buildings")```## Compare centroids in a map### Interactive map```{r mapcomparecentr}#| fig-format: pngmapview(Centroid_geo, col.region = "blue") + mapview(Centroid_pop, col.region = "red") + mapview(Centroid_build, col.region = "black") + mapview(Municipalities_geo, alpha.regions = 0) # polygon limits```See how the building, population and geometric centroids of Sintra are apart, from closer to Tagus, to the rural area.### Static mapTo produce the same map, using only `plot()` and `st_geometry()`, we need to make sure that the geometries have the same crs.quarto-executable-code-5450563D```r#| eval: falsest_crs(Centroid_geo) # 4326 WGS84st_crs(Centroid_pop) # 3763 Portugal TM06```So, we need to transform the geometries to the same crs.quarto-executable-code-5450563D```rCentroid_pop =st_transform(Centroid_pop, crs =4326)Centroid_build =st_transform(Centroid_build, crs =4326)```Now, to use `plot()` we incrementally add layers to the plot.quarto-executable-code-5450563D```r#| fig.dpi: 300#| fig-cap: "Static map of different centroids of Municipalities"#| fig-format: png# Plot the Municipalities_geo polygons first (with no fill)plot(st_geometry(Municipalities_geo), col =NA, border ="black")# Add the Centroids_geo points in blueplot(st_geometry(Centroid_geo), col ="blue", pch =16, add =TRUE) # add!# Add the Centroid_pop points in redplot(st_geometry(Centroid_pop), col ="red", pch =16, add =TRUE)# Add the Centroid_build points in blackplot(st_geometry(Centroid_build), col ="black", pch =16, add =TRUE)```In the [next section](desire-lines.qmd) we will use these centroids to calculate the desire lines between them.quarto-executable-code-5450563D```r#| include: false#| eval: falsest_write(Centroid_pop, "../data/Centroid_pop.gpkg", delete_dsn =TRUE)```
INE. 2022. “Censos 2021- XVI Recenseamento Geral da População. VI Recenseamento Geral da Habitação.” Lisboa: Instituto National de Estatística. https://censos.ine.pt/xurl/pub/65586079.